/*
    Copyright 2008-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Contains the baes_grain_model, which models the grain model used
    by the Baes et al galaxy validation test case. */

// $Id$

#ifndef __baes_grain_model__
#define __baes_grain_model__

#include <vector>
#include "grain_collection_grain_model.h"
#include "grain_size.h"
#include "thermal_equilibrium_grain.h"
#include "misc.h"

namespace mcrx {
  template <template<class> class, typename> class baes_grain_model;
  template <template<class> class chromatic_policy, typename T_rng_policy> 
  class baes_grain_model_creator_functor;
};


/** This class implements the grain model used by the Baes et al
    galaxy validation test case. The class methods are not reentrant
    and can not be used by multiple threads. For the temperature
    calculation, you must make a copy of the object for each
    thread. Note that the baes_grain.fits is an average cross section
    curve, so it doesn't really have a size distribution or numbers of
    grains. The mass opacity normalization is arbitrary (set to be the
    same as the p04 grain), so the amount of dust MUST be set by
    testing the opacity and multiplying to get the desired optical
    depth.
  */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy = mcrx::global_random> 
class mcrx::baes_grain_model : 
  public mcrx::grain_collection_grain_model<chromatic_policy, T_rng_policy> {

  void retarded (const std::string& data_directory, const Preferences& p,
		 const T_unit_map& units);

public:
  baes_grain_model (const std::string&, const Preferences& p, 
		   const T_unit_map&);
  baes_grain_model (const baes_grain_model& g) :
    // need to make sure array members are copied. We make the copies
    // thread-local, since the methods are not reentrant anyway.
    grain_collection_grain_model<chromatic_policy, T_rng_policy>(g) {};

  virtual boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> > 
  clone() const {
    return boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> >
      (new baes_grain_model(*this));};
};


/** This is constructor functionality but is in a distinct function
    since gdb can't debug constructors. */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
void
mcrx::baes_grain_model<chromatic_policy, T_rng_policy>::
retarded (const std::string& data_directory, const Preferences& p,
	  const T_unit_map& units)
{
  // shortcuts
  std::vector<boost::shared_ptr<emitting_grain> >& grains = 
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::grains_;
  std::vector<array_1>& dnda =
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::dnda_;

  // copy units from those supplied into *scatterer* unit map.
  T_unit_map& u = 
    mcrx::scatterer<chromatic_policy, T_rng_policy>::units_;
  u["length"] = units.get("length");
  u["wavelength"] = units.get("wavelength");
  u["mass"] = units.get("mass");

  // The grains use their preferred units, which are also copied into
  // the grain model units.

  // Read the grains.
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + "/baes_grain.fits", p)));

  T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  u2["length"]=grains.front()->units().get("length");
  u2["wavelength"]=grains.front()->units().get("wavelength");
  // grains have no mass unit. But the size distributions do, so for
  // simplicity we select the scatterer mass unit as our mass unit.
  // NOTE that this arrangement means the size distributions and the
  // grains in general use different length units.
  u2["mass"] = u.get("mass");

  // we hardcode "size distribution" here, as it is just the number of
  // grains per mass.

  // this is the number of grains in 1Msun, ~7.6e46
  T_float npermass = constants::Msun/
    (4./3.*constants::pi*0.12e-6*0.12e-6*0.12e-6*3.6e3);
  // apply mass unit
  npermass *= units::convert(units.get("mass"),"Msun");

  // delta_size() of the grain will return 1. The "size distribution"
  // should be in scatterer units, so dnda here should be number of
  // grains per unit mass DIVIDED by the length unit conversion. This
  // will ensure that dn/da*da is n per unit mass.
  
  // evaluate the size distributions. Note that we need to convert
  // from grain length unit here to the scatterer units
  const T_float gsize_conv = units::convert(u2.get("length"),
					    u.get("length"));

  dnda.resize(grains.size());
  dnda[0].resize(grains[0]->asizes().size());
  dnda[0] = npermass/gsize_conv;


}

template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
mcrx::baes_grain_model<chromatic_policy, T_rng_policy>::
baes_grain_model (const std::string& data_directory, const Preferences& p,
		 const T_unit_map& units) :
  grain_collection_grain_model<chromatic_policy, T_rng_policy>()
{
  retarded(data_directory, p, units);
}


template <template<class> class chromatic_policy, typename T_rng_policy> 
class mcrx::baes_grain_model_creator_functor
{
public:
  typedef boost::shared_ptr<mcrx::grain_model<chromatic_policy, T_rng_policy> > return_type;
  
  baes_grain_model_creator_functor() {};
  return_type operator()(const Preferences& p) const {
    const std::string sigma_dir = 
      word_expand(p.getValue("grain_data_directory", std::string ()))[0];
    
    //THIS IS A HACK BECAUSE WE CAN ONLY HAVE ONE ARGUMENT
    T_unit_map units;
    units["mass"] = p.getValue("massunit", std::string());
    units["length"] = p.getValue("lengthunit", std::string());
    units["wavelength"] = p.getValue("wavelengthunit", std::string());

    return return_type
      ( new mcrx::baes_grain_model<chromatic_policy, T_rng_policy> 
	(sigma_dir, p, units));
  };
};


#endif
